clear;
clc;

%% This script generates the figure illustrating how RMSE varies with Tn, nobs, and qn.
% To create the figure:
% - Run main_Server.m for all combinations of the following settings.
% - Save the resulting output files in the ./Results/ folder.

Tnlist = 5:5:50; 
nobslist = 21*[78, 26, 13, 26/3, 13/3, 13/7, 1, 1/3, 1/7];
qnlist = [21, 42, 63];

rmse = nan(length(Tnlist),length(nobslist),length(qnlist));
rmse_tr = nan(length(Tnlist),length(nobslist),length(qnlist));
bias= nan(length(Tnlist),length(nobslist),length(qnlist));
bias_tr= nan(length(Tnlist),length(nobslist),length(qnlist));
std= nan(length(Tnlist),length(nobslist),length(qnlist));
std_tr= nan(length(Tnlist),length(nobslist),length(qnlist));

rmseJ = nan(length(Tnlist),length(nobslist),length(qnlist));
rmse_trJ = nan(length(Tnlist),length(nobslist),length(qnlist));
biasJ= nan(length(Tnlist),length(nobslist),length(qnlist));
bias_trJ= nan(length(Tnlist),length(nobslist),length(qnlist));
stdJ= nan(length(Tnlist),length(nobslist),length(qnlist));
std_trJ= nan(length(Tnlist),length(nobslist),length(qnlist));

for i = 1:length(Tnlist)
disp(Tnlist(i))
    for j = 1:length(nobslist)

        for k = 1:length(qnlist)

           y = load(['./Results/Simu_data_Tn_',num2str(Tnlist(i)),'_nobs_',num2str(nobslist(j)),'_qn_',num2str(qnlist(k)),'.mat']);


           bias(i,j,k) = real(mean(y.Lambda_hat_AJX(1,:) - y.LambdaT_C)); 
           std(i,j,k) = sqrt(var(abs(y.Lambda_hat_AJX(1,:)  - y.LambdaT_C))); 
           bias_tr(i,j,k) = real(mean(y.Lambda_hat_AJX_truebeta(1,:)  - y.LambdaT_C)); 
           std_tr(i,j,k) = sqrt(var(abs(y.Lambda_hat_AJX_truebeta(1,:)  - y.LambdaT_C))); 
           rmse(i,j,k) = sqrt(mean(y.Lambda_hat_AJX(1,:)  - y.LambdaT_C).^2+var(y.Lambda_hat_AJX(1,:)  - y.LambdaT_C));
           rmse_tr(i,j,k) = sqrt(mean(y.Lambda_hat_AJX_truebeta(1,:)  - y.LambdaT_C).^2+var(y.Lambda_hat_AJX_truebeta(1,:)  - y.LambdaT_C));

           biasJ(i,j,k) = real(mean(y.Lambda_hat_AJX(2,:) - y.LambdaT_J)); 
           stdJ(i,j,k) = sqrt(var(abs(y.Lambda_hat_AJX(2,:)  - y.LambdaT_J))); 
           bias_trJ(i,j,k) = real(mean(y.Lambda_hat_AJX_truebeta(2,:)  - y.LambdaT_J)); 
           std_trJ(i,j,k) = sqrt(var(abs(y.Lambda_hat_AJX_truebeta(2,:)  - y.LambdaT_J))); 
           rmseJ(i,j,k) = sqrt(mean(y.Lambda_hat_AJX(2,:)  - y.LambdaT_J).^2+var(y.Lambda_hat_AJX(2,:)  - y.LambdaT_J));
           rmse_trJ(i,j,k) = sqrt(mean(y.Lambda_hat_AJX_truebeta(2,:)  - y.LambdaT_J).^2+var(y.Lambda_hat_AJX_truebeta(2,:)  - y.LambdaT_J));

        end
        

    end


end

%%

figure;
heatmap(Tnlist, {'5 min','15 min','30 min','45 min','90 min','3.5 hr','1 day','3 days'},(round(rmse(:,1:end-1,1),3))','ColorScaling','log')
title('Continuous Risk Premia')
ax = gcf;
exportgraphics(ax,'rmse_C.pdf','Resolution',300) 
%%

figure;
heatmap(Tnlist, {'5 min','15 min','30 min','45 min','90 min','3.5 hr','1 day','3 days'},(round(rmseJ(:,1:end-1,1),3))','ColorScaling','log')
title('Jump Risk Premia')
ax = gcf;
exportgraphics(ax,'rmse_J.pdf','Resolution',300) 